home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_4.3 / xcc / circ_geo.c < prev    next >
C/C++ Source or Header  |  1999-09-11  |  9KB  |  357 lines
  1. /* 
  2.  * circ_geo.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * CIRC(le)_GEOM(etry)
  11.  *
  12.  * collection of routines to evaluate geometrical properties 
  13.  * of circular shapes
  14.  */
  15.  
  16. #include "xcc.h"
  17.  
  18. #define    SYMM        0
  19. #define    ASYMM        1
  20. #define    CCC_MODE    ASYMM          /* mode of center coord computation */
  21.  
  22. /*
  23.  * evaluate slope of linear segment
  24.  */
  25. double
  26. slopev (struct spoint pt1, struct spoint pt2)
  27. {
  28.   double m;
  29.   double max_slope = 500.0;     /* max. poss. slope */
  30.   double r1, c1, r2, c2;
  31.  
  32.   r1 = (double) pt1.y;
  33.   c1 = (double) pt1.x;
  34.   r2 = (double) pt2.y;
  35.   c2 = (double) pt2.x;
  36.  
  37.   if (F_TO_INT (r2 - r1) == 0) {
  38.     if (F_TO_INT (c2 - c1) == 0)
  39.       m = 0.0;
  40.     else
  41.       m = max_slope;
  42.   }
  43.   else
  44.     m = (c2 - c1) / (r2 - r1);
  45.  
  46.   return (m);
  47. }
  48.  
  49.  
  50. /*
  51.  * evaluate moments of order 0 (m00) and of order 1 (m10, m01)
  52.  * 
  53.  * note: not very useful for evaluation of center coord, because disks
  54.  *    are not sampled ``symmetrically'': --> center pos biased
  55.  */
  56. struct spoint
  57. eval_centroid (int *x, int *y, int n)
  58. {
  59.   double xi, ximm, yi, yimm;
  60.   double m00, m00_sum;
  61.   double m10, m10_sum, m01, m01_sum;
  62.   struct spoint ctr;
  63.   int i;
  64.  
  65. /*
  66.  * compute moments
  67.  * (reference: Hu, IRE Trans. Inf. Theory, IT-8, 179-187 (1962) )
  68.  * (see also module pmom.c)
  69.  */
  70.   m00 = m10 = m01 = 0.0;
  71.   for (i = 1; i < n; i++) {
  72.     ximm = (double) (*(x + i - 1));
  73.     xi = (double) (*(x + i));
  74.     yimm = (double) (*(y + i - 1));
  75.     yi = (double) (*(y + i));
  76.  
  77.     m00_sum = 0.5 * (yi * ximm - xi * yimm);
  78.     m00 += m00_sum;
  79.  
  80.     m10_sum = 0.5 * (xi + ximm) * m00_sum;
  81.     m10_sum -= 0.5 * ((yi - yimm) * (xi * xi + xi * ximm + ximm * ximm) / 6.0);
  82.     m10 += m10_sum;
  83.  
  84.     m01_sum = 0.5 * (yi + yimm) * m00_sum;
  85.     m01_sum += 0.5 * ((xi - ximm) * (yi * yi + yi * yimm + yimm * yimm) / 6.0);
  86.     m01 += m01_sum;
  87.   }
  88.   ctr.y = (short) F_TO_INT (m10 / m00);
  89.   ctr.x = (short) F_TO_INT (m01 / m00);
  90.  
  91.  
  92. #ifdef DBG_EVAL_CENTROID
  93.   printf ("EVAL_CENTROID -- centroid of %3d-vertex polygon:", n);
  94.   printf (" (%3d, %3d)\n", ctr.r, ctr.c);
  95. #endif
  96.  
  97.   return (ctr);
  98. }
  99.  
  100.  
  101.  
  102. /*
  103.  * estimate center coords of circular shape based on (horizontal) 
  104.  * disk chord and two line segments by constructing the intersection 
  105.  * of the respective perpendicular bisectors
  106.  *
  107.  * symmetric case (CCC_MODE = SYMM):
  108.  *    construct line segments connecting respective right and left end 
  109.  *    points of disk chord and edge tuple;
  110.  * asymmetric case (CCC_MODE = ASYMM):
  111.  *   construct line segment connecting respective left end points of
  112.  *   disk chord and edge tuple, as well as line segment connecting
  113.  *   left end point of disk chord with right end point of edge tuple
  114.  *   (preferred on basis of error analysis)
  115.  */
  116. double
  117. pbi_ctr (int ir, struct edge_tuple *cetpl, struct triple *cdsk_lchrd)
  118. {
  119.   struct spoint eL, eR;
  120.   struct spoint cL, cR;
  121.  
  122.   double dcc;
  123.   double ecL_r, ecR_r;
  124.  
  125.  
  126. /* initialize variables */
  127.  
  128.   cL.y = cR.y = cdsk_lchrd->r;
  129.   cL.x = cdsk_lchrd->cl;
  130.   cR.x = cdsk_lchrd->cr;
  131.  
  132.   eL.y = eR.y = ir;
  133.   eL.x = cetpl->cl;
  134.   eR.x = cetpl->cr;
  135.  
  136. /*      dcc = (double)cdsk->ctr.x; */
  137.   dcc = 0.5 * (0.5 * (cL.x + cR.x) + 0.5 * (eL.x + eR.x));
  138.  
  139.  
  140.   /* ctr r-coord derived from left edge points */
  141.   ecL_r = 0.5 * (eL.y + cL.y) - (dcc - 0.5 * (eL.x + cL.x)) * slopev (cL, eL);
  142.  
  143.   /* ctr r-coord derived from right edge points */
  144.   if (CCC_MODE == SYMM)
  145.     ecR_r = 0.5 * (eR.y + cR.y) - (dcc - 0.5 * (eR.x + cR.x)) * slopev (cR, eR);
  146.  
  147.   /* ctr r-coord derived from right eTpl and left dChrd edge points */
  148.   if (CCC_MODE == ASYMM)
  149.     ecR_r = 0.5 * (eR.y + cL.y) - (dcc - 0.5 * (eR.x + cL.x)) * slopev (cL, eR);
  150.  
  151.   return (0.5 * (ecL_r + ecR_r));
  152. }
  153.  
  154.  
  155. /*
  156.  * for future use (?) -- to be tested!!
  157.  *
  158.  * construct estimate of center coords on the basis of osculating circles
  159.  * computed from successive triples of vertices
  160.  *
  161.  * for explicit formula, see e.g.:
  162.  * ref: Pavlidis, ``Alg for Graphics and Img Proc'', sect 10, probl 10.11
  163.  *               Comp Sci Press, Rockville, Md, 1982
  164.  */
  165. struct spoint
  166. eval_ctr (struct spoint *v, int n)
  167. {
  168.   struct spoint ctr;
  169.   struct spoint *pvmm, *pv, *pvpp;
  170.  
  171.   double dmr, dpr, dmc, dpc;
  172.   double smr, spr, smc, spc;
  173.   double c1, c2, d;
  174.  
  175.   double cr, scr, cc, scc;
  176.   long nv;
  177.   int i;
  178.  
  179.   scr = scc = 0.0;
  180.   nv = 0;
  181.   for (i = 0; i < n; i++) {     /* form all triples of non-coll vertices */
  182.     pv = &v[i];
  183.     if (i == 0)
  184.       pvmm = &v[n - 1];
  185.     else
  186.       pvmm = &v[i - 1];
  187.     if (i == n - 1)
  188.       pvpp = &v[0];
  189.     else
  190.       pvpp = &v[i + 1];
  191.  
  192.     dmr = pvmm->y - pv->y;      /* x1-x2 */
  193.     dmc = pvmm->x - pv->x;      /* y1-y2 */
  194.     dpr = pv->y - pvpp->y;      /* x2-x3 */
  195.     dpc = pv->x - pvpp->x;      /* y2-y3 */
  196.  
  197.     if (fabs (d = (dmr * dpc - dpr * dmc)) < 0.000001) {
  198.       d = 0.0;
  199.       printf ("EVAL_CTR -- WARNING: vertices collinear\n");
  200.     }
  201.     else {
  202.       smr = pvmm->y + pv->y;    /* x1+x2 */
  203.       smc = pvmm->x + pv->x;    /* y1+y2 */
  204.       spr = pv->y + pvpp->y;    /* x2+x3 */
  205.       spc = pv->x + pvpp->x;    /* y2+y3 */
  206.  
  207.       c1 = dmr * smr + dmc * smc;
  208.       c2 = dpr + spr + dpc * spc;
  209.  
  210. /* evaluate sum of ctr coords, to be averaged */
  211.       cr = 0.5 * (c1 * dpc - c2 * dmc) / d;
  212.       scr += cr;
  213.       cc = 0.5 * (c2 * dmr - c1 * dpr) / d;
  214.       scc += cc;
  215.       nv++;
  216.  
  217.       printf ("EVAL_CTR -- ");
  218.       printf ("nv=%ld, cr=%lf, cc=%lf\n", nv, cr, cc);
  219.     }
  220.   }
  221.   ctr.y = (short) F_TO_INT (scr / nv);
  222.   ctr.x = (short) F_TO_INT (scc / nv);
  223.  
  224.  
  225. #ifdef DBG_EVAL_CTR
  226.   printf ("EVAL_CTR -- center of circ, based on %3d-vertex polygon:", n);
  227.   printf (" (%3d, %3d)\n", ctr.r, ctr.c);
  228. #endif
  229.  
  230.   return (ctr);
  231. }
  232.  
  233.  
  234. /*
  235.  * find best circle to fit set of points {(x_i, y_i), 0<=i<n}
  236.  * with center position given by centroid: compute radius of inertia
  237.  *
  238.  * ref: Pavlidis, ``Alg for Graphics and Img Proc'', sect 12.6.2, probl 12.5
  239.  *               Comp Sci Press, Rockville, Md, 1982
  240.  */
  241. int
  242. oeval_rad (struct spoint *v, int n, struct disk *cdsk)
  243. {
  244.   double sdrsq, sdcsq;
  245.   double r_bar, c_bar;
  246.   int i;
  247.  
  248.   r_bar = (double) cdsk->ctr.y;
  249.   c_bar = (double) cdsk->ctr.x;
  250.   sdrsq = sdcsq = 0.0;
  251.   for (i = 0; i < n; i++) {
  252.     sdrsq += SQR ((v + i)->y - r_bar);
  253.     sdcsq += SQR ((v + i)->x - c_bar);
  254.   }
  255.  
  256.   return (F_TO_INT (sqrt ((sdrsq + sdcsq) / n)));
  257. }
  258.  
  259.  
  260. /*
  261.  * find best circle to fit set of points {(x_i, y_i), 0<=i<n}
  262.  * with center position given by centroid: compute radius of inertia
  263.  *
  264.  * ref: Pavlidis, ``Alg for Graphics and Img Proc'', sect 12.6.2, probl 12.5
  265.  *               Comp Sci Press, Rockville, Md, 1982
  266.  */
  267. int
  268. eval_rad (struct disk *cdsk)
  269. {
  270.   struct spoint *v;
  271.   double sdrsq, sdcsq;
  272.   double r_bar, c_bar;
  273.   int ich, nch;
  274.   int iv, nv;
  275.  
  276.   nch = cdsk->nch;
  277.   if ((v = (struct spoint *) calloc (2 * (nch + 1), sizeof (struct spoint))) == NULL) {
  278.     printf ("\n...mem alloc for array of vertices v failed!!\n");
  279.     return (0);
  280.   }
  281.  
  282. /*
  283.  * extract vertex coords
  284.  */
  285.   for (ich = 0; ich < nch; ich++) {
  286.     (v + ich)->y = (v + ((2 * nch - 1) - ich))->y = cdsk->chords[ich].r;
  287.  
  288.     (v + ich)->x = cdsk->chords[ich].cl;
  289.     (v + ((2 * nch - 1) - ich))->x = cdsk->chords[ich].cr;
  290.   }
  291.  
  292.  
  293.   r_bar = (double) cdsk->ctr.y;
  294.   c_bar = (double) cdsk->ctr.x;
  295.   sdrsq = sdcsq = 0.0;
  296.   nv = 2 * nch;
  297.   for (iv = 0; iv < nv; iv++) {
  298.     sdrsq += SQR ((v + iv)->y - r_bar);
  299.     sdcsq += SQR ((v + iv)->x - c_bar);
  300.   }
  301.   free (v);
  302.  
  303.   return (F_TO_INT (sqrt ((sdrsq + sdcsq) / nv)));
  304. }
  305.  
  306.  
  307.  
  308. /*
  309.  * generate a new estimate of the best circle for the current set
  310.  * of contour points: estimate the center in form of the centroid,
  311.  * then compute a (suboptimal) value of the radius in form of the
  312.  * radius of inertia
  313.  */
  314. int
  315. oeval_circ (struct disk *cdsk)
  316. {
  317.   struct spoint *v;
  318.   int ich, nch;
  319.  
  320.   nch = cdsk->nch;
  321.   if ((v = (struct spoint *) calloc (2 * (nch + 1), sizeof (struct spoint))) == NULL) {
  322.     printf ("\n...mem alloc for array of vertices v failed!!\n");
  323.     return (0);
  324.   }
  325.  
  326. /*
  327.  * evaluate new centroid position and update entry in disk structure
  328.  */
  329.   for (ich = 0; ich < nch; ich++) {
  330.     (v + ich)->y = (v + ((2 * nch - 1) - ich))->y = cdsk->chords[ich].r;
  331.  
  332.     (v + ich)->x = cdsk->chords[ich].cl;
  333.     (v + ((2 * nch - 1) - ich))->x = cdsk->chords[ich].cr;
  334.   }
  335.  
  336. #ifdef DBG_EVAL_CTR
  337.   printf ("EVAL_CIRC -- polygon vertex coords, {(r, c}\n");
  338.   for (ich = 0; ich < 2 * nch; ich++)
  339.     printf ("  (%3d, %3d)\n", (v + ich)->y, (v + ich)->x);
  340. #endif
  341.  
  342.   cdsk->ctr = eval_ctr (v, 2 * nch);
  343.   printf ("\n...new estimated center position: (%3d, %3d)\n",
  344.           cdsk->ctr.y, cdsk->ctr.x);
  345.  
  346. /*
  347.  * estimate new radius
  348.  */
  349.   cdsk->rad = oeval_rad (v, 2 * nch, cdsk);
  350.   printf ("...estimated radius: %d\n", cdsk->rad);
  351.  
  352.  
  353.   free (v);
  354.  
  355.   return (1);
  356. }
  357.